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P^ Abstract The relativistic motion of a compact binary system moving in circular 

^_^ orbit is investigated using the post-Newtonian (PN) approximation and the pertur- 

^ bative self-force (SF) formalism. A particular gauge-invariant observable quantity 

'""' is computed as a function of the binary's orbital frequency. The conservative ef- 

VO feet induced by the gravitational SF is obtained numerically with high precision, 

^~~^ and compared to the PN prediction developed to high order. The PN calculation 

i__i involves the computation of the 3PN regularized metric at the location of the par- 

O tide. Its divergent self-field is regularized by means of dimensional regularization. 

">-<* The poles °^ (d — 3)^^ which occur within dimensional regularization at the 3PN 

^^ order disappear from the final gauge-invariant result. The leading 4PN and next-to- 

leading 5PN conservative logarithmic contributions originating from gravitational- 
wave tails are also obtained. Making use of these exact PN results, some previously 
unknown PN coefficients are measured up to the very high 7PN order by fitting to 
the numerical self-force data. Using just the 2PN and new logarithmic terms, the 
value of the 3PN coefficient is also confirmed numerically with very high precision. 
\^ The consistency of this cross-cultural comparison provides a crucial test of the very 

psj different regularization methods used in both SF and PN formalisms, and illustrates 

J. I the complementarity of these approximation schemes when modelling compact bi- 

^^ nary systems. 
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1 Introduction and motivation 

For the gravitational wave observatories LIGOA'^irgo/GEO on Earth and LISA 
in space, routine identification of inspiralHng compact binaries (binary systems 
composed of neutron stars and/or black holes) will require high-accuracy predic- 
tions from general relativity theory |1]J|2]. Providing such predictions represents a 
formidable task that can be addressed using approximation schemes in general rel- 
ativity. The two main approximation schemes available are: (i) the post-Newtonian 
expansion, well suited to describe the inspiralling phase of compact binaries in the 
slow motion and weak field regime independently of the mass ratio, and (ii) the self- 
force approach, based on perturbation theory, which gives an accurate description 
of extreme mass ratio binaries even in the strong field regime. 

The post-Newtonian (PN) templates for compact binary inspiral have been devel- 
oped to 3.5PN order in the phase |[3l|4]|3 and 3PN order in the amplitude IHEl (see 
Blanc het's contribution in this volume)|^ These are suitable for the inspiral of two 
neutron stars in the frequency bandwidth of LIGO and Virgo detectors. For detec- 
tion of black hole binaries (with higher masses), the construction of template banks 
either requires the matching of the PN waveform with full numerical simulations 
for the merger phase and the ringdown of the final black hole 1,8, ,9J, or using the 
effective-one-body formalism ifTOl (see also Nagar's contribution in this volume). 

In a completely different parameter regime, gravitational self-force (SF) analysis 
ifTTl [T2I [131 [T4I \T5i (see also Poisson's contribution in this volume) is expected 
to provide templates for extreme mass ratio inspirals (EMRIs) anticipated to be 
present in the LISA frequency bandwidth. SF analysis is a natural extension of first 
order perturbation theory, and the latter has a long history of comparisons with post- 
Newtonian analysis QglElIISlIlllEOlEllEllEl. SF analysis itself, however, is 
just now mature enough to present some limited comparisons with PN analysis, 
although it is not yet ready for template generation. 

In recent works ll24l I25JI (hereinafter referred to as Papers I and II respectively) 
we performed a high-accuracy comparison between the PN and SF analyses in their 
common domain of validity, that of the slow motion weak field regime of an ex- 
treme mass ratio binary (see illustration of various methods in Fig. [Til. The problem 
was tackled previously by Detweiler fSSl, who computed numerically within the 
SF approach a certain gauge invariant quantity (called the redshift observable), and 
compared it with the 2PN prediction extracted from existing PN results ll27l . We 
then extended this comparison in Papers I-II up to higher PN orders. This required 
an improvement in the numerical resolution of the SF calculation in order to dis- 
tinguish more accurately the various contributions of very high PN order terms. 
However, our primary difficulty has been that the relevant PN results for the metric 
were previously not available beyond the 2PN order, and had to be carefully derived. 
We have finally demonstrated an excellent agreement between the SF contribution 



^ As usual the nPN order refers to terms equivalent to (v/c)^" beyond Newtonian theory, where v 
is a typical internal velocity of the material system and c is the speed of light. 
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to the analytical PN result (derived through 3PN order, with inclusion of specific 
logarithmic terms at 4PN and 5PN orders) and the exact numerical SF result. 
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Fig. 1 Different analytical approximation schemes and numerical techniques are used to study 
black hole binaries, depending on the mass ratio mi /mi and the orbital velocity i'^ ~ Gm/ri2, 
where m = mi +;7!2- The post-Newtonian theory and black hole perturbation theory can be com- 
pared in the slow motion regime (v<^ f equivalent tori 2 2> Gm/c~ for circular orbits) of an extreme 
mass ratio (mi <Si mi) binary. 



In this article we present a summary of the Papers I and II. The plan is as follows: 
After having introduced in Sec. |2] the coordinate-invariant relation used to perform 
the comparison, we describe in Seep] how the divergent self-field of point parti- 
cles is regularized in both SF and PN formalisms. Sec.fflprovides a brief overview 
of how the SF computation proceeds. Sees. l5] and l6] present the post-Newtonian 
computations of the 3PN regularized metric, as well as of the 4PN and 5PN loga- 
rithmic contributions. The PN results for the gauge-invariant relation are discussed 
in Sec. 17] Finally, Sec. [8] is devoted to the comparison of these PN results with the 
SF numerical data, and the measurement of unknown high order PN coefficients. 



4 Luc Blanchet, Steven Detweiler, Alexandre Le Tiec and Bernard F. Whiting 

2 The gauge invariant redshift observable 

We consider a system of two (non-spinning) compact objects with masses mi and 
m2, and moving on slowly inspiralling quasi-circular orbits. In the PN analysis, let 
mi and m2 be arbitrary; in the SF analysis, further assume that nii ^ m2- We can 
then call nii the "particle", and m2 the "black hole". 

Self-force analysis shows that the dissipative parts of the self-force for a circu- 
lar orbit are the t and (p components. These result in a loss of energy and angular 
momentum from the small mass at the same precise rate as energy and angular 
momentum are radiated away [26 1. In addition, earlier perturbative calculations of 
energy and angular- momentum fluxes lHH [iTl HH HH llO] H] |22l |23l for fliis situ- 
ation show them to be equivalent to the results of the PN analysis in their common 
domain of validity. Hence, by invoking an argument of energy and angular momen- 
tum balance, we know that the PN results also agree with the dissipative parts of the 
SF, and further comparison can reveal nothing new. 

For our PN-SF comparison, we shall thus neglect the dissipative, radiation- 
reaction force responsible for the inspiral, and restrict ourselves to the conservative 
part of the dynamics. In PN theory this means neglecting the dissipative radiation- 
reaction force at 2.5PN and 3.5PN orders, and considering only the conservative 
dynamics at the even-parity IPN, 2PN and 3PN orders. This clean separation be- 
tween conservative even-parity and dissipative odd-parity PN terms is correct up to 
3.5PN order However such split breaks at 4PN order, since at that approximation 
arises a contribution of the radiation-reaction force, which originates from grav- 
itational wave tails propagating to infinity {2K\ (this will be further discussed in 
Sec.|6]l. In SF theory there is also a clean split between the dissipative and conser- 
vative parts of the self-force (see e.g. ||29]| ). This split is particularly transparent for 
a quasi-circular orbit, where the r component is the only non-vanishing component 
of the conservative self- force. 

Henceforth, the orbits of both masses are assumed to be and to remain circular, 
because we are ignoring the dissipative radiation-reaction effects. For our compari- 
son we require two physical quantities which are precisely defined in the context of 
each of our approximation schemes. The orbital frequency Q of the circular orbit as 
measured by a distant observer is one such quantity. The second quantity is defined 
as follows. 

With circular orbits and no dissipation, the geometry has a heUcal Killing vector 
field A:". A Killing vector is only defined up to an overall constant factor In our case 
A:" extends out to a large distance where the geometry is essentially flat. There, 

k"da = d, + Qd^, (I) 



in any natural coordinate system which respects the helical symmetry llSOll . We let 
this equality define the overall constant factor, thereby specifying the Killing vector 
field uniquely. 

An observer moving with the particle mi, while orbiting the black hole m2, would 
detect no change in the local geometry. Thus the four-velocity u" of the particle 
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is tangent to the Killing vector k" evaluated at the location of the particle, which 
we denote by k". A second physical quantity is then defined as the constant of 
proportionality, call it uf , between these two vectors, namely 

u'^^u^kf. (2) 

The four-velocity of the particle is normalized so that {ga[i)]'^i'^\ = ~1; {8ap)] is 
the regularized metric at the particle's location, whereas the metric itself is formally 
singular at the particle mi in both PN and SF approaches. 

If we happen to choose a coordinate system such that ([T]) is satisfied everywhere, 
then in particular k[~ I, and thus uj = u\, the t component of the four- velocity of 
m\. The Killing vector on the particle is then k" — u" /u[, and simply reduces to the 
particle's ordinary post-Newtonian coordinate velocity vf /c. In such a coordinate 
system, the description of the invariant quantity we are thus considering is 



aJ\-l/2 



V'l V'l 



(3) 



It is important to note that this quantity is precisely defined in both PN and SF 
frameworks, and it does not depend upon the choice of coordinates or upon the 
choice of perturbative gauge. The quantity u^ represents the redshift of light rays 
emitted from the particle and received on the helical symmetry axis perpendicular 
to the orbital plane [26 1; we shall refer to it as the redshift observable. 



3 Regularization issues in the SF and PN formalisms 

The redshift observable (|3]l depends upon using a valid method of regularization. 
The regularized metric {g^p ) , is defined with very different prescriptions in the SF 
and PN approaches. Both analyses require subtle treatment of singular fields at the 
location of the masses. Subtracting away the infinite part of a field while carefully 
preserving the part which is desired is always a delicate task. 
In the SF prescription, the regularized metric reads 

8api^)^8ap{x)+hlpi^)^ (4) 

where g^p denotes the background Schwarzschild metric of the black hole, and 
where the "Regular" metric perturbation h^n is smooth in a neighborhood of the 
particle, and given by the difference 

'^lp='Cp~l4ip (5) 

between the retarded metric perturbation /z^l and the purely locally determined 
"Singular" field h^g . Following the DetweUer- Whiting prescription lfT3l . a Hadamard 
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expansion of Green's functions in curved spacetime provides an expansion for h^„ 
lfT3l . In a neighborhood of the particle with a special, locally-inertial coordinate 
system, h^„ appears as the mi/r part of the particle's Schwarzschild metric along 
with its tidal distortion caused by the background geometry of the large black hole. 
Details of the expansion are given in Sec. 6.1 of OTI . The special locally inertial 
coordinates for a circular geodesic in the Schwarzschild metric are given as func- 
tions of the Schwarzschild coordinates in Appendix B of |32|. (See also Detweiler's 
and Barack's contributions in this volume.) Since the metric (|4]) is regular at the 
particle's position y", we simply have 

i8fp),-8%iyi)- (6) 

In the perturbative SF analysis we are only working through first order in the mass 
ratio q = mi/m2, and at that level of approximation /i^g = i^{q). Then uj can be 
computed accurately to the same perturbative order and compares well with the post- 
Newtonian result to 2PN order ll26l . The regularized 2PN metric is known |27 1, and 
therefore the comparison is straightforward. 

In the PN prescription on the other hand, one first computes the metric §|j^g(x,f) 
using an iterative PN procedure at any field point outside the particle, in a coordi- 
nate system x" — {ct,x'}. That metric is generated by the two particles, and includes 
both regular and singular contributions around each particle. Such iterative PN cal- 
culation is a very long and intricate procedure up to say 3PN, at which order it 
will be partly based on existing computations of the 3PN equations of motion using 
Hadamard |33| and dimensional |34| regularizations. Then we compute the regular- 
ized metric at the location of the particle by taking the limit when x — ?► yi (f ), where 
yi(?) is the particle's trajectory. In 3 spatial dimensions, that limit is singular In 
order to treat the infinite part of the field, we extend the computation in d spatial di- 
mensions, following the prescription of dimensional regularization 1351 [36l . which 
is based on an analytic continuation (AC) in the dimension d viewed as a complex 
number. We do not use Hadamard's regularization which found its limit at 3PN or- 
der. Considering the analytic continuation in a neighborhood of £ = t/ — 3 — > 0, we 
define 



(^™„) =AC lim/„^«(x,f) 



1 £^0 



,PN, 



x^yi 



(7) 



The limit e ^ does not exist in general due to the presence of poles oc e^' occuring 
at the 3PN order; we thus do not take the strict limit e ^- but compute the singular 
Laurent expansion when e ^- 0. At 3PN order the result takes the schematic form 

(/a^)i-^4'^(ynO+4(yi,0 + ^(e), (8) 

where g^a (yi , t) denotes the pole part, which is purely of 3PN order, and g^i (yi , f ) 
is the finite part. At higher PN orders we expect the presence of multipole poles 
oc £^". As we shall see, the (simple) poles at 3PN order will disappear from the 
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final gauge invariant relationship uJ{Q). In fact the occurence of poles at the 3PN 
order is specific to the use of the harmonic gauge condition. Previous work on the 
3PN equations of motion of point particle binaries has shown that the poles can be 
absorbed into a renormalisation of the worldlines of the particles, so they should 
not appear in any physical coordinate invariant quantity. Thus dimensional regu- 
larization is a powerful regularization method in the PN context. In particular this 
regularization is free of the ambiguities plaguing the Hadamard regularization at the 
3PN order Il37ll34l l5l (see the contributions by Schafer and Blanchet in this volume). 
Although the two regularizations in SF and PN analyses have been carefully de- 
signed, it appears to be non trivial that they will yield results consistent up to a high 
level of approximation. Our cross cultural comparison of the redshift observable uj 
is a test of the equivalence of the SF and PN metrics (|6| and ^ and is, thus, a test 
of the two independent (and very different) regularization procedures in use. 



4 Circular orbits in the perturbed Schwarzschild geometry 

Previously we described the truly coordinate and perturbative-gauge independent 
properties of Q and the redshift observable u^ . In this section we use Schwarzschild 
coordinates, and we refer to "gauge invariance" as a property which holds within 
the restricted class of gauges for which ([T]l is a helical Killing vector. In all other 
respects, the gauge choice is arbitrary. With this assumption, no generality is lost, 
and a great deal of simplicity is gained. 

The effect of the gravitational self-force is most easily described as having mi 
move along a geodesic of the regularized metric (|4]). We are interested in circular 
orbits and let m" be the four-velocity of miplThis differs from the four-velocity u" 
of a geodesic of the straight Schwarzschild geometry at the same radial coordinate 
r by an amount of ^(^). Recall that we are describing perturbation analysis with 
q-^1, therefore /i^g — &'{q), and all equations in this Section necessarily hold only 
through first order in q. 

It is straightforward to determine the components of the geodesic equation for 
the metric Q [26], and then to find the components of the four-velocity m" of m\ 
when it is in a circular orbit at Schwarzschild radius r. We reiterate that the four- 
velocity is to be normalized with respect to g^n +/^a6 rather than ^„n, and that /z^g 
is assumed to respect the symmetry of the helical Killing vector In this case we 
hav£] 



{u'f = -^—- [l +M«m'5/,R _ 'a^uf^drhl. 



r-3m2L- «P 2 ^ ^^ 



(«") 



2 



r 



2m2 pW2(l+5"M^/i^^) 1_„_^ R 

> , V It U Or't' rj f- 



r{r-3m2) 



r{r — 2m2) 2 



(9a) 
(9b) 



Since we are interested in the motion of the small particle m i , we remove the index 1 from u " . 
' In all of this section we shall set G = c = 1 . 
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A consequence of these relations is that the orbital frequency of mi in a circular 
orbit about a perturbed Schwarzschild black hole of mass OT2 is, through first order 
in the perturbation, given by 

The angular frequency 12 is a physical observable, and is independent of the gauge 
choice. However the perturbed Schwarzschild metric does not have spherical sym- 
metry, and the radius of the orbit r is not an observable and does depend upon the 
gauge choice. That is to say, an infinitesimal coordinate transformation of ff{q) 
might change u"uPdrh^g. But if it does, then it will also change the radius r of the 

orbit in just such a way that Q^ as determined from ( [TO| i remains unchanged. Both 
u' = u^ and u'^ = Qu^ are gauge invariant as well. 

Our principle interest is in the relationship between Q and u^, which we now 
establish directly using ( |9a] l and ( [TOJ I. To this end, following ll26l . we introduce the 
gauge invariant measure of the orbital radius 

and readily establish a first order, gauge invariant, algebraic relationship between u^ 
(to which u' evaluates in our gauge) and Rq (or equivalently Q), namely: 

(u^f={i-'^y\i+a"u^hlp). (12) 

See Paper I for a detailed derivation of this result. The lowest order term in q on 
the right-hand-side is identical to what is obtained for a circular geodesic of the 
unperturbed Schwarzschild metric. Indeed, recall that the Schwarzschild part of u^ 
is known exactly as u^^^^ = (1 ^ ^'ni/^^n)^ ■ Thus, if we write 

«^ = «Schw+^«SF + ^(?^): (13) 

the first order term in ( fTS] ) gives: 

?4f=^1-^)" •i"'i^h%, (14) 

which is ^{q), and contains the effect of the "gravitational self-force" on the re- 
lationship between u^ and Q, even though it bears little resemblance to a force. 
The numerical SF approach henceforth focuses attention on the calculation of the 
combination m"m^/i^„. See |l26l| and Paper I for more details on the implementation 
of the regularization of the perturbation, and on the numerical computation of the 
quantity u"uPh^„. 
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5 Overview of the 3PN calculation 

Our aim is to compute the 3PN regularized metric (|7| by direct post-Newtonian 
iteration of the Einstein field equations in the case of singular point mass sources. 
In the dimensional regularization scheme, we look for the solution of the Einstein 
field equations in d+l space-time dimensions. We treat the space dimension as 
an arbitrary complex number, li e C, and interpret any intermediate formula in the 
PN iteration of those equations by analytic continuation in d. Then we analytically 
continue d down to the value of interest (namely 3), posing d = 3 + e. In most of 
the calculations we neglect terms of order £ or higher, i.e. we retain the finite part 
and the eventual poles. 



5.1 Iterative PN computation of the metric 

Defining the gravitational field variable h"^ = \/—gg"^ — 77"",rjand adopting the 
harmonic coordinate condition dnh"l^ — 0, we can write the "relaxed" Einstein field 
equations in the form of ordinary d'Alembert equations, namely 

Uh^P ^l^^lglT'^P +A"P[h,dh,d^h], (15) 

where D = r\^^ dj^dy is the^af-spacetime d'Alembertian operator in d+\ dimen- 
sions. The gravitational source term A "^ in ( [TS] ) is a functional of h^^ and its first 
and second space-time derivatives; it depends explicitly on the dimension d. The 
matter stress-energy tensor T"^ is composed of Dirac delta-functions in d dimen- 
sions. Finally, the li-dimensional gravitational constant G^''' is related to the usual 
Newton constant G by 

Gt'" = G.0 

where (.q denotes the characteristic length associated with dimensional regulariza- 
tion. We shall check that this length scale disappears from the final gauge-invariant 
3-dimensional result. 

The 3PN metric is given in expanded form for general matter sources in terms 
of some "elementary" retarded potentials (sometimes called near-zone potentials) 
y, Vi, K, Wij, Ri, X, Zij, % and f , which were introduced in Ref. ||33| for 3 di- 
mensions and generalized to d dimensions in Ref. |34|. All these potentials have a 
finite non-zero post-Newtonian limit when c -^ +oo and parameterize the successive 
PN approximations. Although this decomposition in terms of near-zone potentials 
is convenient, such potentials have no physical meaning by themselves. 

Let us first define the combination 



g''''=G£o, (16) 



* Here g"P is the contravariant metric, inverse of the covariant metric g^^ of determinant 
g = det(g„^), and 7]°" = diag(— 1, 1, 1, 1) represents an auxiliary Minkowski metric in Cartesian 
coordinates. 
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2 /d-3\ AX 16f 



■ = y- 



d-2 



K 



(17) 



Then the 3PN metric components can be written in the rather compact form 



PN 

Sm 



PN 



„-2r/, 



(d-i)f 

{d-2)c- 




-10 



), 



(18a) 



^PN^ ,rf-2).2 
"7 



4 . 



16 

~6 



2,j-ViVj 



8Ri 16 

c^ c' 


^■ + 


^%-^. 




(18b) 


2(^-2) ^"''^'^ 




c-^), 



(18c) 



where the exponentials are to be expanded to the order required for practical calcula- 
tions. The successive PN truncations of the field equations ( [T5] l give us the equations 
satisfied by all the above potentials up to 3PN order We conveniently define from 
the components of the matter stress-energy tensor T"P the following density, current 
density, and stress density: 



2 (t/-2)r'*' + r" 

d-l c^ 



CTi 



a,j^T'\ 



(19a) 

(19b) 
(19c) 



where T" = dijT'-i . As examples, the leading-order potentials in the metric obey 



8, J 



(^kk 



d-2 



1 U-l 



2\d-2 



diVdjV . 



(20a) 
(20b) 

(20c) 



All the potentials evidently include many PN corrections. The potentials V and V, 
have a compact support (i.e. their source is localized on the isolated matter sys- 
tem) and will admit a finite limit when £ — > without any pole. Most of the other 
potentials have, in addition to a compact-support part, a non-compact support con- 
tribution, such as that generated by the term oc diVdjV in the source of Wij. These 
non-compact support pieces are the most delicate to compute, because they typically 
generate some poles oc 1/e at the 3PN order The d'Alembert equations satisfied by 
all higher-order PN potentials, whose sources are made of non-linear combinations 
of lower-order potentials, can be found in Paper 1. Clearly, the higher the PN order 
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of a potential, the more complicated is its source, but it requires computations at a 
lower relative order 

Many of the latter potentials have already been computed for compact binary 
systems, and we have extensively used these results from If33ll34l . Notably, all the 
compact-support potentials such as V and V,-, and all the compact-support parts of 
other potentials, have been computed for any field point x, and then at the source 
point yi following the regularization. However, the most difficult non-compact sup- 
port potentials such as X and T could not be computed at any field point x, and 
were regularized directly on the particle's world-line. Since for the equations of 
motion we needed only the gradients of these potentials, only the gradients were 
regularized on the particle, yielding the results for (<3,^)i and {diT)^ needed in the 
equations of motion. However the 3PN metric requires the values of the potentials 
themselves regularized on the particles, i.e. (X) j and (T) j . For the present work we 
therefore computed, using the tools developed in ll33l[34l . the difficult non-linear 
potentials {X) j and (T) j , and especially the non-compact support parts therein. Un- 
fortunately, the potential X is always the most tricky to compute, because its source 
involves the cubically-non-linear and non-compact-support term Wij djjV, and it has 
to be evaluated at relative IPN order. 

In this calculation we also met a new difficulty with respect to the computation 
of the 3PN equations of motion. Indeed, we found that the potential X is divergent 
because of the bound of the Poisson-like integral at infinity. Thus the potential X de- 
velops an IR divergence, in addition to the UV divergence due to the singular nature 
of the source and which is cured by dimensional regularization. The IR divergence 
is a particular case of the well-known divergence of Poisson integrals in the PN 
expansion for general (regular) sources, linked to the fact that the PN expansion is 
a singular perturbation expansion, with coefficients typically blowing up at spatial 
infinity. The IR divergence is discussed in Paper I, where we show how to resolve it 
by means of a finite part prescription. 

The 3PN metric ([T8]l is valid for a general isolated matter system, and we apply 
it to the case of a system of A^ point-particles with "Schwarzschild" masses m„ and 
without spins (here a= 1 , • • • , A^). In this case we have 

C7(x,r)=£A«5('')[x-y„(f)], (21a) 

a 

cT,(x,o=E^"^«^*''H'^-y«(0]> (21b) 

a 

a,j{x, t)=Y^ii, vi vi, 5^'''^ [x - y„(f )] , (21c) 

a 

where d^'^' denotes the Dirac density distribution in d spatial dimensions, such that 
/d''x5''''(x) = 1. We defined the effective masses of the particles by 

M„(f) = , "'" == , (22a) 
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2 



i"n(0 = 



d-\ 



d-2- 



V21 



Ai«(0- (22b) 



5.2 The example of the zeroth-order iteration 

For illustration purposes, let us consider the simple case of the dimensional reg- 
ularization of the Newtonian potential generated by two point particles of masses 
m\ and m2- The PN metric ( fTSJ l reduces to ^™ = — 1 +2V /c^ + ^(c^^), where 
y is solution of the Poisson equation AV = —'XnG^'^^G + i^{c^^), with source 

(7(x) = ^^J^[mi5('')(x-yi)+'«25('')(x-y2)] +^(c-2). Solving this Poisson 
equation in d space dimensions yields 

2{d-2) / G^^^'mi G^''hn2 \ , ^, _2^ 

"^^^ V|x-yiK-2 + R^y^ 



y(x).^^^^M^:V^ + ^^;£^)+^(c-^), (23) 



where A: = F(^^)/7rT" tends to 1 when £ — > (F is the usual Eulerian function). 
When d — 3, the Newtonian potential ( |23| ) is not defined in the limit x — > yi be- 
cause of the divergent self-field of particle 1 . By contrast, thanks to the analytic 
continuation in the space dimension, it is always possible to choose 3i{d) < 2 such 
that the li-dimensional potential ( |23l ) has a well-defined limit when x — > yi, namely 
{V)i = y(yi) given by 



2{d-2) ,_G(^)m2 

yd-'. 

12 



{V),^^j-^k-^ + ^{c-'). (24) 



where we have posed ri2 = |yi — y2|- Relying on the unicity of the analytic contin- 
uation, we obtain the unique 3 -dimensional result 

(V)i = — + ^(e) + ^(c-2). (25) 

'•12 

This procedure is clear up to a high order, but let us mention a subtle point in 
the calculation of Paper I, namely that we had to systematically re-introduce the 
correction terms ^(e) in the Newtonian part of the metric and other quantities. 
Indeed, in various operations such as replacing ri2 in Eq. ( p5| l by its 3PN expression 
in terms of the orbital frequency Q, the corrections €?(£) will be multiplied by some 
poles at 3PN order, and they will therefore contribute in fine to the finite part at 3PN 
order Such corrections are necessary only in the Newtonian results (since the poles 
arise only at 3PN order). For instance, the 3-dimensional Newtonian potential ( |25] l 
is to be replaced by its li-dimensional version valid up to terms i^(e^), namely 



ri2 



2 V 4 



+ &{e^) \ + ^ic-"-) , (26) 
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in which p = \/~^e^'^ with C = 0.5772 • • • being the Euler-Mascheroni constant. 
This Newtonian calculation is generalized at higher orders in the iterative PN 
process described in the previous Section. The result is the 3PN regularized metric 
([8]l which is given in explicit form by Eqs. (4.2) of Paper I. 



6 Logarithmic terms at 4PN and 5PN orders 

We now discuss the logarithmic contributions in the near-zone metric of a generic 
isolated post-Newtonian source, and then of a compact binary system. Our motiva- 
tion is that knowing analytically determined logarithmic terms in the PN expansion 
is crucial for efficiently extracting from the SF data the numerical values of higher 
order PN coefficients. This will be further discussed in Sec.|8] 

The occurence of logarithmic terms in the PN expansion has been investigated 
in many previous works |l38l|39l|40l|4T]|42l|43 44, 28, 45|. Notably Anderson et 
al. BTl found that the dominant logarithm arises at the 4PN order, and Blanchet & 
Damour f^S", ^45 1 showed that this logarithm is associated with gravitational wave 
tails modifying the usual 2.5PN radiation-reaction damping at the 4PN order. Fur- 
thermore the general structure of the PN expansion is known [1441 : it is of the type 
£(v/c)*^[ln(v/c)]''', where k and q are positive integers, involving only powers of 
logarithms; more exotic terms such as [ln(ln(v/c))]'' cannot arise. 

Following Paper II, we shall determine the leading 4PN logarithm and the next- 
to-leading 5PN logarithm in the conservative part of the dynamics of a compact 
binary system. The computation of such logarithmic contributions relies only very 
weakly on a regularization scheme. We shall thus work in 3 space dimensions; from 
now on we set e = 0. 



6.1 Physical origin of logarithmic terms 

Because of the non-linearity of the field equations, the gravitational field at coordi- 
nate time t is in general not a function of the state of motion of the source at retarded 
time t — r/c, where r — |x| is the distance to the center of the source, but depends 
on the entire past "history" of the source. This means that the near-zone metric de- 
pends on the source at all times before the current time f, say t' ^t (indeed in the 
near-zone one expands all retardations, i.e. r/c -^ 0). This "hereditary" effect starts 
at 4PN order in the near-zone, and originates from gravitational-wave tails, namely 
the scattering of gravitational radiation by the background curvature of the space- 
time generated by the mass M of the source. In a certain gauge where the hereditary 
terms are collected into the 00 component of the metric, we have 

«'(».') - -'g?.."..71-/«P'('')ln (^) , (27, 
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(n) 

where M is the ADM mass of the source, and M^^ is the n-th time derivative of its 
quadrupole moment. This tail-induced contribution is of 4PN order. 

The occurence of the tail effect in the near-zone metric implies that the usual 
2.5PN radiation-reaction force density in the matter source is corrected at the rela- 
tive 1 .5PN order as 



7.C 



(5), , 4GM f . n-), „ f cit 



21 



(28) 



where p is the Newtonian mass density in the source, and X is the typical wavelength 
of the radiation. The leading term in ( p8] l is the standard Burke-Thorne f^H" "471 
radiation-reaction potential at the 2.5PN order, responsible for the leading radiation 
effect. The hereditary correction was obtained in Il28ll45]| . and shown to be consis- 
tent with wave tails propagating at large distances from the source f^Sl. 

The radiation-reaction force (|28| deserves its name because it is not invariant 
under a time reversal, and therefore gives rise to dissipative effects. A good way to 
see this is to change the condition of retarded potentials to advanced potentials, i.e. 
to formally change c into — c. The first term in ( |28| l is clearly non invariant because it 
comes with an odd number of powers of 1/c. The second term is also non invariant, 
despite the fact that it comes with an even power of 1/c in front (i.e. 1/c^ in ( [28] ), 
corresponding to 4PN), because it is composed of an integral extending over the 
past rather than a time-symmetric integral. 

However, thanks to the even power of 1/c carried by the hereditary integral pT] ), 
this means that there exists a conservative piece associated with it. Recall that in our 
calculation we neglect the dissipative radiation-reaction effects and are interested 
only in the conservative part of the dynamics; we have implemented this restriction 
by assuming the existence of the helical Killing vector Q, depending on the orbital 
frequency Q of the circular motion. The presence of this scale Q , imposed by the 
helical Killing symmetry, permits immediately to identify the conservative piece 
associated with the tail term in ( |27] i, through the decomposition 

K^)-"(^)-'"(i)' 

where now the characteristic length scale X is defined by A = 2kc/Q. The second 
term in (|29|, when inserted into the tail integral in ( pTJ i, can be integrated out and 
gives rise to the conservative 4PN piece 

5^or(x,0 = ^-^^"^'Mf^{t) In (£) . (30) 

This contribution to the metric has to be included in our study of the conservative 
part of the dynamics, while the purely dissipative piece in ( |27| i enters the radiation 
reaction ( |28] l, and is excluded by our assumption of existence of the helical Killing 
vector. The term pO] ) is precisely the conservative 4PN logarithmic contribution in 
the near-zone metric that we want to compute, as well as its 5PN correction. 
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Notice that it is possible to find a gauge where all the logarithms are gathered in 
the 00 component of the metric, as we did in Eq. ( |30l ), only at 4PN order When look- 
ing to subdominant logarithmic terms at 5PN order, we are obliged to include some 
vectorial 0/ and tensorial ij components of the metric. But still it will be possible, 
and extremely convenient, to define a gauge in which the logarithms are "maxi- 
mally" transferred to the 00 and 0/ components (while the remaining contribution 
in the ij components is "minimal"). Using such a gauge saves a lot of calculations 
when performing the PN iteration; this was the strategy adopted in Paper II to com- 
pute the higher order 5PN logarithms beyond Eq. (l30l. 



6.2 Expression of the near-zone metric 

To compute these 4PN and 5PN conservative logarithmic contributions, we make 
use of the multipolar post-Minkowskian wave-generation formalism |44 28l l48l 
I45JI . We first identify these logarithmic contributions in the exterior of a generic 
isolated post-Newtonian source. We then deduce the metric inside the matter source 
by a matching performed in the exterior part of the near-zone. We finally get the 
4PN and 5PN logarithmic contributions in the near-zone metric, valid in a specific 
gauge defined in Paper II, asj 



5^00 = 



^10 



1-: 



2U 



x-yf 



35c2 



..vyf 






5 ci2 "* 



^gij 



G^-M 

^10 



d^x' 



x-x' 



p'x'^'ln 



■ e 



16 

21" 



i-y^ 



^^ ^"ab "'J 



-e .x"^5<^' 



In^ 



In 



(i) 



(i) 



,12 



,14 



189c2 






In 



(D 

(31a) 
(31b) 
(31c) 



where M^h, Mahc and Sab denote the mass quadrupole, mass octupole and current 
quadrupole moments of the source respectively, and where 



f/(x,f)==G 



d^x' 



x-x' 



^P(x',0 



(32) 



is the Newtonian potential, sourced by the Newtonian mass density p of the source. 
We recall that A = 2kc/Q, where Q. is the scale entering the Killing vector ([T]). 
Notice in particular the 5PN contribution in 5goo which involves the Poisson integral 
of a logarithmically modified source density, and which will contribute in fine to the 
5PN logarithms in the case of binary systems. 



^ We use shorthands such as x"* = x"x^ ; x"'"-' = . 



(S" 



■ S'^x* -I- 5'"-Y')r^ denotes the 



symmetric and trace-free part of x" ; e^bc is the Levi-Civita antisymmetric symbol. 
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We can then apply the result ( (3T] i to the case of a compact binary system moving 
on a circular orbit. In this case, Q is the orbital frequency of the motion. We compute 
the metric at the location of one of the particles and obtain the result 



(5goo)i = 





[8/ 2f/r 

[5 Iv' c2 ^ 


)ytMfrl, 


u^yh'iM't^ 




1 ^ v2v"^Af(**' ^ ^flbc.A^) 

35c2^^^^ "'' 189c2^i "'"■ 


■"(T)-<i) 


, (33a) 


G'-M 


'Y/Ml'-%e.byrs^S 


■"(T)-<i)- 


(33b) 


G'-M 


'\yfMf^^, 


ln(^).. 


( 1 
^Cl2 


/ 


(33c) 



(5goOi = 



(%;) 



where L'^i = Gm2/ri2 is the Newtonian potential felt by particle 1. Finally, the last 
step is to replace the multipole moments Mab, Mate and Sab by the relevant PN ex- 
pressions valid for circular-orbit compact binaries; we need also the ADM mass M, 
which reduces to m = mi +m2 in first approximation. Note that the mass quadrupole 
moment Mab (and also the ADM mass M) must crucially include a IPN contribu- 
tion. Again, we emphasize that for the 4PN and 5PN logarithms we do not need 
the full apparatus of dimensional regularization, in contrast to the fully fledged 3PN 
calculation sketched in Sec.|5] 



7 Post-Newtonian results for the redshift observable 



To compute the gauge invariant quantity u^ (associated with the particle 1 for helical 
symmetry, circular orbits), we adopt its coordinate form as given by ([3]l, namely 



?a(3j 



v'^vl 



-1/2 



(34) 



and plug into it the 3PN regularized metric (|8]l, explicitly computed from the 3PN 
near-zone expression ([TSj reduced to binary point masses, and including the 4PN 
and 5PN logarithmic corrections computed in Sec. [6] To begin with, this yields 
the expression of u' for an arbitrary mass ratio q = nii/mi, and for a generic non- 
circular orbit in a general reference frame. We then choose the frame of the center of 
mass, which is consistently defined by the nullity of the center-of-mass integral of 
the motion, deduced from the equations of motion. Restricting ourselves to exactly 
circular orbits (consistently with the helical Killing symmetry we neglect radiation- 
reaction effects), the result is expressed by means of the convenient dimensionless 
gauge invariant PN parameter 
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(35) 



which is directly related to the orbital frequency £1 of the circular orbit, and depends 
on the total mass m = mi+m2of the binary. 

We discover most satisfactorily that all the poles °^ l/e (as well as the associated 
constant £q) cancel out in the final expression for u'. Our final result for a 3PN (plus 
4PN and 5PN logarithmic terms), gauge invariant, algebraic relationship between 
u^ (to which u' now evaluates) and x (or equivalently Q), ian 
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(36) 



We introduced the notation A = (7«2 — m\)/m = ^/\ — 4v, where v = 11111112/ m^ is 
the symmetric mass ratio. WTiile it has been shown in ll26l (see also Sec. [2] above) 
that M^ is gauge invariant at any PN order in the extreme mass ratio limit V ^ 1, here 
we find that it is also gauge invariant for any mass ratio up to 3PN order (even up to 
5PN order for the logarithmic terms). This result is expected from (|2]), according to 
which u^ is a scalar under our hypothesis of helical symmetry. Being proportional 
to the symmetric mass ratio v, the 4PN and 5PN logarithmic contributions vanish in 
the test-mass limit — this is clear given that the Schwarzschild result for u^{Q) does 
not involve any logarithm. Notice that the functions A4(v) and As{v) entering the 
expression of the non-logarithmic contribution to u^{Q) at the 4PN and 5PN orders 
are unknown, and would be very difficult to compute within standard PN theory. 
However we know that they are polynomials in v, with leading-order coefficient 
given by the Schwarzschildean result [see Eqs. (|40ji]. 

We now investigate the small mass ratio regime ^ <C 1, for comparison purposes 
with the perturbative SF calculation. We introduce a convenient PN parameter ap- 
propriate to the small mass limit of particle 1 ; 



' The Landau o symbol for remainders takes its standard meaning. 
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which is related to the usual PN parameter xhy x = y{l + q)^'^, and to the gauge- 
invariant measure ( [TT| ) of the orbital radius hy y = Gm2/{RQC^). We also use the 
expression of the symmetric mass ratio v in terms of the (asymmetric) mass ratio 
q — m\/m2, namely v = q/{l +q)^. Our complete redshift observable, expanded 
through post-self-force order, is of the type 

"^ = "Schw + 1 «SF + 1^ «PSF + ^{q^), (38) 

where the Schwarzschild result is known in closed form as Mgchw ~ (^ ^ -^y) . By 
expanding the PN result ( |36) l in powers of q, we find that the self-force contribution 
reads 

4f(3') = -y - 2/ - 5/ + (^- 1^ + 11 7^2^ / 

a4 - y \ny) / + (as + ^ \ny) / + o(/) . (39) 

The coefficients a4 and a^ are pure numbers which parametrize the small mass ratio 
expansions of the functions A4 and A5 through 

15309 / 25515\ ^, 2, 

168399 / 168399 \ ^, 2^ ,,^,., 

We also give the result for the combination u^uf^h^n related to u^p by Eq. ( [I4| , 
since this is the quantity primarily used in the numerical SF calculation: 

a BfR . 273/ 1387 41 2\ 4 

u^uPhlp = -2y-y^ _ _/ + (^___ + _^2j ^4 

fl4 - IP Iny) / + («5 + ^ Iny) / + o(/) . (41) 

We have conveniently rescaled the first-order perturbation h^„ by the mass ratio q, 

denoting h^o = /i^g /q. Here 04 and as denote some unknown pure numbers related 
to a^ and a^ by 

9301 123 2 

fl4 = 2a4 + — — -— -7r% (42a) 

64 il 



17097 369 



,2 



«5 = 2a5-3a4 + ^^-— ;r^ (42b) 
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The expansions (|39l)-(|4T]) were determined up to 2PN order °^ y^ in ll26l . based on 
the Hadamard-regularized 2PN metric given in ll27l . The result at 3PN order °^ y^ 
was obtained in Paper I using the powerful dimensional regularization scheme. By 
comparing the expansion ( (39] l with our accurate numerical SF data for u^p{Q), we 
shall be able to measure the coefficients a4 and aj (or 04 and 05) with at least 8 
significant digits for the 4PN coefficient, and 5 significant digits for the 5PN coeffi- 
cient. These results, as well as the estimation of even higher-order PN coefficients, 
will be detailed in the next Section. 

Similarly, from the PN result (|36]l valid for any mass ratio q, we get the post-self- 
force contribution as 

r / N o 2 97 3 /725 41 jA 4 

+ e4y'+U + ^lny\y'^ + oiy'), (43) 

which could in principle be compared to a future post-self-force calculation making 
use of second-order black hole perturbation theory. Note that there is no logarithm 
at 4PN order in the post-self-force term; the next 4PN logarithm would arise at 
cubic order q^, i.e. at the post-post-SF level. The coefficients £4 and £5 in ( |43] l are 
unknown, and unfortunately they are expected to be extremely difficult to obtain, 
not only analytically in the standard PN theory, but also numerically as they require 
a second-order perturbative SF scheme. 



8 Numerical evaluation of post- Newtonian coefficients 

In the self-force limit, the SF effect u^p on the redshift observable u^ is related via 
J 12b to the regularized metric perturbation h^o at the location of the particle through 

4p=i(l-3y)-i/2 «««-%, (44) 

where S" is the background four- velocity of the particle. Recall that here /i^„ stands 
for the perturbation per unit mass ratio, that is h^a/q- In SF analysis, the combi- 
nation U^Uni^n arises more naturally than Mgp; this is the quantity we shall be in- 
terested in fitting in this Section. However our final results in Table l5] will include 
the corresponding values of the coefficients for the redshift observable u^^. We refer 
to Sec. II of Paper I for a discussion of the computation of the regularized metric 
perturbation h^n, and the invariant properties of the combination M^u^/i^n with re- 
spect to the choice of perturbative gauge. In this Section we often use r = 1/y, a 
gauge invariant measure of the orbital radius scaled by the black hole mass m2 [see 
Eqs. ([Tl)and([37]l]. 
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Our earlier numerical work in ll26l l24l I25JI provided values of the function 
u"uni^g{r) which cover a range in r from 4 to 750. Following a procedure de- 
scribed in 1 32 1, we have used Monte Carlo analysis to estimate the accuracy of our 
values for u"uni^„. As was reported in Paper I, this gives us confidence in these 

base numbers to better than one part in 10^^. We denote a standard error a repre- 
senting the numerical error in S^w^/i^g by 

ac:i\u"ul^hlp\xExlO-^\ (45) 

where E ~ 1 is being used as a placeholder to identify our estimate of the errors in 
our numerical results. 



8.1 Overview 

A common task in physics is creating a functional model for a set of data. In our 
problem we have a set of A^ data points /, and associated uncertainties a,, with 
each pair evaluated at an abscissa r,. We wish to represent this data as some model 
function f{r) which consists of a linear sum of M basis functions Fj{r) such that 

M 

f{r) = Y.^jFj('-)- (46) 

.7=1 

The numerical goal is to determine the M coefficients Cj which yield the best fit in a 
least squares sense over the range of data. That is, the Cj are to be chosen such that 



N 

(=1 



fi-LUcjF.in) 



(47) 



is a minimum under small changes in the Cj. For our application we choose the 
basis functions Fj{r) to be a set of terms which are typical in PN expansions, such 
as r^', r^^, . . . , and also terms such as r^^ lii(r). We recognize that a solution to 
this extremum problem is not guaranteed to provide an accurate representation of 
the data (r,,/,, a,). The quality of the numerical fit is measured by x^ as defined in 



Eq. ( |47| . If the model of the data is a good one, then the x^ statistic itself has an 
expectation value of the number of degrees of freedom in the problem, N —M, with 
an uncertainty (standard deviation) of \J2{N — M). 

Our numerical work leans heavily upon Ref. |49J for solving the extremum prob- 



lem for Eq. (47 1. The numerical evaluation of the fitting coefficient Cj includes a 
determination of its uncertainty T,j which depends upon i) the actual values of r, in 
use, ii) all of the a,, and iii) the set of basis functions Fj{r). In fact, the estimates of 
the I.J do not depend at all on the data (or residuals) being fitted. As a consequence 
the estimates of the I,j are only valid if the data are well represented by the set of 
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basis functions. For emphasis: the Xj depend upon Fj{ri) and upon O; but are com- 
pletely independent of the ft. Only if the fit is considered to be good, could the T,j 
give any kind of realistic estimate for the uncertainty in the coefficients Cj. If the 
fit is not of high quality (unacceptable x^), then the Ej bear no useful information 
B9l . We will come back to this point in the discussion below. 

We also should remark that the task of determining coefficients in the 1 /r charac- 
terization of our numerical data is almost incompatible with the task of determining 
an asymptotic expansion of u"'uPh^g from an analytic analysis. Analytically, the 
strict r — > +00 limit is always technically possible, whereas numerically, not only is 
that limit never attainable, but we must always contend with function evaluations 
at just a finite number of discrete points, obtained within a finite range of the inde- 
pendent variable, and computed with finite numerical precision. Nevertheless, this 
is what we have done below. 

The numerical problem is even more constrained than we have just indicated. At 
large r, even though the data may still be computable there, the higher order terms 
for which we are interested in evaluating PN coefficients rapidly descend below 
the error level of our numerical data. This is clearly evident in Fig. |2] below. For 
small r, the introduction of so many PN coefficients is necessary that it becomes 
extremely difficult to characterize our numerical data accurately. Thus, in practice, 
we find ourselves actually working with less than the full range of our available 
data. At large r we could effectively drop points because they contribute so little to 
any fit we consider At the other extreme, the advantage of adding more points in 
going to smaller r is rapidly outweighed by the increased uncertainty in every fitted 
coefficient. This results from the need to add more basis functions in an attempt to 
fit the data at small r. Further details will become evident in Sec.l8.4lbelow. 



8.2 Framework for evaluating PN coefficients numerically 

In a generic fashion we describe an expansion of u'^rrh^a in terms of PN coeffi- 
cients Uj and bj with 

«"«-%- L;|fr-in^L;|fr> (48) 

whsTe flo is the Newtonian term, ai is the IPN term and so on. Similarly, for use 
in applications involving m^ we also introduce the coefficients aj and j3y in the 
expansion of the SF contribution 

4p=E|^-ln.E|^. (49) 

These series allow for the possibility of logarithmic terms, which are known not to 
start before the 4PN order. We also concluded in Paper II that (Inr)^ terms cannot 
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arise before the 5.5PN order. Since we are computing a conservative effect, possible 
time-odd logarithmic squared contributions at the 5.5PN or 6.5PN orders do not 
contribute. But there is still the possibility for a conservative VPN (Inr)^ effect, 
probably originating from a tail modification of the dissipative 5.5PN (Inr)^ term. 
However, we shall not permit for such a small effect in our fits. As discussed below 
in Sec. |8.4[ we already have problems distinguishing the VPN linear Inr term from 
the VPN non-logarithmic contribution. 

The analytically determined values of the coefficients aq, ai, 02, 03, b4, bs and 
Oo, «!, a2, 0:3, j34, jSs computed in Ref. [26 1 and Papers I-II are reported in Tablefl] 



coeff. 


value 


coeff. 


value 


"0 


_2 


Oo 


-1 
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-1 
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-2 


<32 
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a2 
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1387 , 41 2 
24 ^ \6^ 


a-i 


_12i , 41 2 
1 + i2^ 


b4 
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Pi 


_64 
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, 5944 


ft 


,^6 



Table 1 The analytically determined PN coefficients for u^uPfi^g (left) and Mgp (right). 



8.3 Consistency between analytically and numerically determined 
PN coefficients 



In this Section we investigate the use of our data for u"ul^h^„ and the fitting proce- 
dures we have described above (and expanded upon in the beginning of Sec. 8.4 1 



We will begin by fitting for enough of the other PN coefficients to be able to verify 
numerically the various coefficients aj,, Ba and b^ now known from PN analysis. 

As a first step in this Section, we will complete the task that was begun in ||26]| . 
namely, the numerical determination of the coefficient 03 (and a^), this time taking 
fully into account the known logarithmic terms at 4PN and 5PN order For illus- 
trative purposes only, these results are given in Table l2] We were able to obtain a 
fit with six undetermined parameters, and could include data from r — VOO down 
to r = 35. Note that, with the inclusion of the b^ and bs coefficients, the precision 
of our tabulated value for 03 has increased by more than four orders of magnitude 
from Paper I, although our accuracy is still no better than about 2Z. Such a discrep- 
ancy is not uncommon. The uncertainty, Z, reflects only how well the data in the 
given, finite range can be represented by a combination of the basis functions. It is 
not a measure of the quality of a coefficient when considered as a PN expansion 
parameter, which necessarily involves an r — > +00 limiting process. 

Our next step is to include the known value for (33 and to use our numerical data 
to estimate values for the b^ and ^75 coefficients. Our best quality numerical result 
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3PN coeff. 


Ref. |26| 


Paper I 


Paper II 


PN (exact) 


a-i 


-32.34(6) 
-27.61(3) 


-32.479(10) 
-27.677(5) 


-32.5008069(7) 
-27.6879035(4) 


-32.50080538- •• 
-27.68790269- •• 



Table 2 The results of a numerical fit for a set of six coefficients that includes a-i , which is now 
known analytically] 24 1 . This fit uses the known results for bn and i>5 1 25 1 , but not the known value 
of a^. Thus, it is not the best-fit of our data possible. The uncertainty in the last digit or two is in 



parentheses. The range runs from r = 
264. 



35 to r = 700, with 266 data points and a respectable x of 



was obtained with five fitted parameters, over a range from r = 700 down to only 
r — 65, and is given in the first row of Table l3] Notice that while our b^ is determined 
relatively precisely, it has only about 6T. accuracy. The higher order coefficient fog 
is more difficult to obtain and, at this point, it is very poorly determined, but we can 
use the known value of bA, in order to improve the accuracy for bs- These results 
are presented in Table |3] which again shows that we needed to fit for a total of six 
parameters to get a result of reasonable accuracy. With this, we have reached a limit 
for treating our data in this way, since adding further parameters and inner points 
does not result in any higher quality fit. 



'"min 


;t-/dof 


04 


b. 


"5 


h 


(26 


bb 


«7 


65 
85 
65 
40 


0.961 
0.976 
0.961 
0.969 


-121.40(1) 
-121.3180(7) 
-121.313(1) 
-121.3052(6) 


-25.612(2) 


-102(1) 
-91.5(7) 
-79(2) 
-47(1) 


45.5(3) 
48.5(2) 
50.6(4) 
55.7(2) 


-2081(9) 
-2170(8) 
-1868(44) 
-359(41) 


131(21) 
625(15) 


-7722(162) 



Table 3 The numerically determined PN coefficients for u"u^h^„ . Each row represents a different 
fit. The first two columns give the starting point fmin at the inner boundary of the fitting range, and 
the value of x~ statistic per degree of freedom (dof) for the chosen fit. The degrees of freedom, 
N — M, for the fit, range between 212 and 255, depending on rmin. If a value for a coefficient is not 
shown, then either that parameter was not included in that particular fit (far right) or its analytically 
known value was used (e.g., bi). The formal uncertainty of a coefficient in the last digit or two is 
in parentheses. The outer boundary is at 700 in each case. 



By now we have presented enough to show that we have data which allows high 
precision, with an accuracy that we now have some experience in relating to the 
computed error estimates. This experience will be valuable when we come to dis- 
cuss further results in the next Section. For convenience, we summarize the relevant 
information further, in Table |4j referring just to our estimates of known PN param- 
eters, and relating our error estimates to the observed accuracy. 
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Paper I 
Table 
Table 3 
Table 3 



coeff. 



hi 



estimate 



-27.677(5) 
-32.5008069(7) 

-25.612(2) 
+55.7(2) 



accuracy 



(11) 
(15) 
(12) 
(9) 



exact result 



-27. 6879- • • 
-32.50080538- 
-25.6 
+56.6095--- 



Table 4 Comparing the analytically known PN coefficients (column 5) with their numerically 
determined counterparts (column 3), and comparing the numerically determined eiTor estimates 
(column 3) with the apparent accuracy (column 4). The source of the data is given in column 1. 



8.4 Determining higher order PN terms numerically 

In this Section we make maximum use of the coefficients which are akeady known. 
We find that in our best fit analysis we can use a set of five basis functions corre- 
sponding to the unknown coefficients 04, 05, og, b(, and aj. 



coeff. 


value 


coeff. 


value 


04 


-121.30310(10) 


04 


-114.34747(5) 


"5 


-42.89(2) 


«5 


-245.53(1) 


"6 


-215(4) 


06 


-695(2) 


be 


+680(1) 


ft 


+339.3(5) 


a-i 


-8279(25) 


Ctl 


-5837(16) 



Table 5 The numerically determined values of higher-order PN coefficients for u^ifih^o (left) and 
for Mjp (right). The uncertainty in the last digit or two is in parentheses. The range runs from r = 40 
to r = 700, with 261 data points being fit. The x^ statistic is 259. We believe that a contribution 
from a b-j confounds the a-j coefficient. Both terms fall off rapidly and have influence over the fit 
only at small r. And the radial dependence of these two terms only differ by a factor of In r [or 
possibly (Inr)^] which changes slowly over their limited range of significance. 



In Table [5] we describe the numerical fit of our data over a range in r from 40 
to 700. The x^ statistic is 259 and slightly larger than the degrees of freedom, 256, 
which denotes a good fit. Further, we expect that a good fit would be insensitive to 
changes in the boundaries of the range of data being fit, and we find, indeed, that if 
the outer boundary of the range decreases to 300 then essentially none of the data 
in the Table changes, except for x^ and the degrees of freedom which decrease in 
a consistent fashion. Figure 2 shows that in the outer part of the range u"uPh^n is 
heavily dominated by only a few lower order terms in the PN expansion — those 
above the lower black double-dashed line in the figure. 

The inner edge of the range is more troublesome. The importance of a given 
higher order PN term decreases rapidly with increasing r. Moving the inner bound- 
ary of the range outward might move a currently well determined term into insignif- 
icance. This could actually lead to a smaller x'^, but it would also lead to an increase 
in the Zj of every coefficient. Moving the inner edge of the range inward might re- 
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25 



1e-08 



r''(uuh-5PNL) 
5PN 




6PN 

7PN 

err 

res I 



100 



200 



300 
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700 



Fig. 2 The absolute value of the contributions of the numerically determined post-Newtonian terms 
to r^u"u" h^g. Here PNL refers to just the logarithm term at the specified order. The contribution 

of 04 is not shown but would be a horizontal line (since the 4PN terms behaves like r^^ ) at approxi- 
mately 121.3 . The remainder after 04 and all the known coefficients are removed from r^u^uPh^g 
is the top (red) continuous line. The lower (black) dotted line labelled "err" shows the uncertainty 
in r^u"uPh^g, namely 2Er^ x 10^'-'. The jagged (green) line labelled "|res|" is the absolute re- 
mainder after all of the fitted terms have been removed. The figure reveals that, with regard to the 
uncertainty of the calculated j("m'^/i^„, the choice £ ~ 1 was slightly too large. 



quire that an additional higher order term be added to the fit. This extra term loses 
significance quickly with increasing r so the new coefficient will be poorly deter- 
mined and also result in an overall looser fit with an increase of Ej for all of the 
coefficients. If the inner boundary and the set of basis functions are chosen prop- 
erly, then a robust fit is revealed when the parameters being fit are insensitive to 
modest changes in the boundaries of the range. The fit described in Table [5] appears 
to be robust. The parameters in this Table are consistent with all fits with the inner 
boundary of the range varying from 35 to 45 and the outer boundary varying from 
300 to 700. 

If an additional term, with coefficient bj, is added to the basis functions then, for 
identical ranges, each of the Ej increases by a factor of about ten, and the changes 
in 04 and 05 are within this uncertainty. The coefficient ae changes sign and b^ and 
aj change by an amount significantly larger than the corresponding Ej. And the new 
coefficient bj is quite large. In the context of fitting data to a set of basis functions 
these are recognized symptoms of over-fitting and imply that the extra coefficient 
degrades the fit. 
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8.5 Summary 



Our best fit can be visualized in Fig.BJ where we plot the self- force effect «gp on the 
redshift variable u^ as a function of r ~ y^^, as well as several truncated PN series 
up to 7PN order, based on the analytically determined coefficients summarized in 
Table [T] as well as our best fit of the higher-order PN coefficients reported in Table 
[5] Observe in particular the smooth convergence of the successive PN approxima- 
tions towards the exact SF results. Note, though, that there is still a small separation 
between the VPN curve and the exact data in the very relativistic regime shown at 
the extreme left of Fig. |3] 




Fig. 3 The self-force contribution Mgp to the redshift observable u^ , plotted as a function of the 
gauge invariant variable y. Note that y^^ is an invariant measure of the orbital radius scaled by the 
black hole mass 1112 [see Eqs. ^TTJ and l |37[ l] . The "exact" numerical points are taken from Ref . 1 26 1 . 
Here, PN refers to all terms, including logarithms, up to the specified order (however recall that we 
did not include in our fit a log-term at 7PN order). 



We have found that our data in the Umited range of 35 ^ r ^ 700 can be ex- 
tremely well characterized by a fit with five appropriately chosen (basis) functions. 
That is, the coefficients in Table|5]are well determined, with small uncertainties, and 
small changes in the actual details of the fit result in coefficients lying within their 
error estimates. Fewer coefficients would result in a very poor characterization of 
the same data while more coefficients result in large uncertainties in the estimated 
coefficients, which themselves become overly sensitive to small changes in specific 
details (such as the actual choice of points to be fitted). In practice, over the data 
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range we finally choose, and with the five coefficients we fit for, we end up with ex- 
ceedingly good results for the estimated coefficients, and with residuals which sink 
to the level of our noise. We have a very high quality fit which is quite insensitive to 
minor details. Nevertheless, as TableHlhints, error estimates for these highest order 
coefficients should be regarded with an appropriate degree of caution. 
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